#include <gtk/gtk.h>
#include "matio.h"
#define w 10	//每个patch有w*w个点
#define n 40	//每个IMAGE中有n*n个patch
#define x 5
#define k (x*x)	//把IMAGE的patch分为k类
#define m 5000
#define size 5	//一个像素点的大小
#define classMax (10*n*n/k)
#define ImageWidthStribe (n*n/k/10)
#define ImageHeightStribe 2
#define ImageWidth (x*(n*n/k*size+ImageWidthStribe))
#define ImageHeight (x*(10*size+ImageHeightStribe))

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct Patch
{
	double data[w*w];
	double value;
};

struct Class
{
	int index[classMax];
	double value;
	int num;
};

//临时变量
int i,j,h,a,b,c,d,t,p;
//分配10*512*512*double = 2.5M*double = 20M，存储十幅图片
double ***data;
//分配patch
struct Patch patch[10*n*n];
//分配class
struct Class class[k];
//标志位
int isChange;
double Max;
double Min;
//所有class中最大中间和最小的Value
double maxValue;
double midValue;
double minValue;
//像素位
guchar rgbbuf[ImageWidth * ImageHeight];
//Widget
GtkWidget *window, *darea;
//Expose Event
GdkEvent event;

/**图像数据处理函数**/
void gdk_data_image()
{
	guchar *pos;
	pos = rgbbuf;
	for(i = 0;i < k;i++)
		class[i].num = 0;
	for (i = 0; i < ImageHeight; i++)
    {
		for (j = 0; j < ImageWidth; j++)
		{
			a = i%(ImageHeight/x) - ImageHeightStribe;
			b = i/(ImageHeight/x);
			c = j%(ImageWidth/x) - ImageWidthStribe;
			d = j/(ImageWidth/x);
			if(a >= 0 && c >= 0 && class[b * x + d].num < classMax)
			{
				double s = patch[class[b * x + d].index[class[b * x + d].num]].value;
				*pos++ = (int)((s - Min)/(Max - Min) * 256);
				class[b * x + d].num++;
			}
			else
				*pos++ = 0;
		}
	}
}

/**kmeans函数**/
void kmeans()
{
	int i,j;
	for(i = 0;i < k;i++)
	{
		double s;
		s = 0;
		for(j = 0;j < classMax;j++)
			s += patch[class[i].index[j]].value;
		s = s/classMax;
		class[i].value = s;
	}
}

/**qsort的比较函数**/
int comp(const void *a,const void *b)
{
	double num = *(double *)b - *(double *)a;
	if(num > 0)   return 1;
	if(num < 0)   return -1;
	return 0;	
}

/**kmedoids函数**/
void kmedoids()
{
	int i,j;
	double temp[classMax];
	for(i = 0;i < k;i++)
	{
		for(j = 0;j < classMax;j++)
			temp[j] = patch[class[i].index[j]].value;
		qsort(temp,classMax,sizeof(double),comp);
		class[i].value = temp[classMax/2];
	}
}

gboolean on_darea_expose (GtkWidget *widget, GdkEventExpose *event, gpointer user_data)
{
	gdk_draw_gray_image (widget->window, widget->style->fg_gc[GTK_STATE_NORMAL],
              0, 0, ImageWidth, ImageHeight, GDK_RGB_DITHER_NORMAL, rgbbuf, ImageWidth);
	return TRUE;
}

int main(int argc, char *argv[])
{
	/**绘图初始化**/
	gtk_init (&argc, &argv);
	window = gtk_window_new (GTK_WINDOW_TOPLEVEL);
	darea = gtk_drawing_area_new ();
	gtk_widget_set_size_request (darea, ImageWidth, ImageHeight);
	gtk_container_add (GTK_CONTAINER (window), darea);
	gtk_signal_connect (GTK_OBJECT (darea), "expose-event", GTK_SIGNAL_FUNC (on_darea_expose), NULL);
	gtk_widget_show_all (window);
	/**数据变量**/
	data = malloc(sizeof(double**)*10);	
	for(i = 0;i < 10;i++)
	{
		data[i] = malloc(sizeof(double*)*512);
		for(j = 0;j < 512;j++)
			data[i][j]= malloc(sizeof(double)*512);
	}
	Max = -1;
	Min = 1;
	maxValue = 0;		
	midValue = 0;
	minValue = 0;
	event.type = GDK_EXPOSE;
	/**打开IMAGES.mat文件**/	
	mat_t *mat = Mat_Open("IMAGES.mat", MAT_ACC_RDONLY);
	//mat_t *mat = Mat_Open("IMAGES_RAW.mat", MAT_ACC_RDONLY);
	if (mat == NULL)
	{
		printf("Error! Cannot Open File!\n");
		return(EXIT_FAILURE);
	}
	/**从IMAGES.mat中读取数据，并且保存到data里面**/
	matvar_t *matvar = Mat_VarReadNextInfo(mat);	
	Mat_VarReadDataAll(mat,matvar);
	for(i = 0;i < 10;i++)
	{
		for(j = 0;j < 512;j++)
		{
			for(h = 0;h < 512;h++)
			{
				data[i][j][h] = ((double*)matvar->data)[512*512*i+512*j+h];
				if(data[i][j][h] > 	Max)
					Max = data[i][j][h];
				if(data[i][j][h] <	Min)
					Min = data[i][j][h];
			}
		}
	}
	/**对每个IMAGE数据进行采样,采样结果为w*w的Patch，每个IMAGE有n*n个Patch**/	
	for(i = 0;i < 10;i++)
	{
		for(j = 0;j < 512;j++)
		{
			a = j / (512 / n);//第i幅图中第a行的patch
			b = j % (512 / n);//该patch中第b行的元素
			if(b < w)
			{
				for(h = 0;h < 512;h++)
				{
					c = h / (512 / n);//第i幅图中第a行第c个的patch
					d = h % (512 / n);//该patch中第b行第d个元素
					if(d < w)
					{	
						t = n*n*i+a*n+c;
						if(t < 10*n*n)
						{
							patch[t].data[b*w+d] = data[i][j][h];
						}
					}
				}
			}
		}
	}	
	/**计算每个patch的平均值**/
	for(i = 0;i < 10*n*n;i++)
	{
		double s;
		s = 0;
		for(j = 0;j < w*w;j++)
			s += patch[i].data[j];
		s = s/(w*w);
		patch[i].value = s;
	}
	/**随机取k个patch作为class的初始值**/	
	for(i = 0;i < k;i++)
		class[i].value = patch[random()%(10*n*n)].value;
	p = 0;
	/**kmeans or kmedoids**/	
	do
	{
		p++;
		/**class状态清空**/
		for(i = 0;i < k;i++)
			class[i].num = 0;
		/**每个patch归属于新的class**/
		for(i = 0;i < 10*n*n;i++)
		{
			/**寻找距离最小的class**/
			double minSubValue = 10;
			int index;
			for(h = 0;h < k;h++)
			{
				/**如果该class还没有满**/							
				if(class[h].num < classMax)
				{
					double s;
					s = fabs(patch[i].value - class[h].value);
					if(s < minSubValue)
					{
						minSubValue = s;
						index = h;
					}					
				}
			}
			if(class[index].num < classMax)
			{
				class[index].index[class[index].num] = i;
				class[index].num++;
			}				
		}
		/**更新image**/
		gdk_data_image();
		gtk_widget_send_expose(window,&event);
		/**计算每个class的value**/
		//kmeans();
		kmedoids();
		/**如果value发生变化了，那么标志isChange为1，表示发生过改变**/
		double temp[k],sMax,sMid,sMin;
		for(i = 0;i < k;i++)
			temp[i] = class[i].value;
		qsort(temp,k,sizeof(double),comp);
		sMax = temp[k-1];
		sMid = temp[k/2];
		sMin = temp[0];
		if(sMid != midValue || sMax != maxValue || sMin != minValue)
			isChange = 1;
		else
			isChange = 0;
		/**把本次的value保存下来，以便下次迭代比较时使用**/
		//if(p % 100 == 0)
			printf("%d %lf %lf %lf\n",p,minValue,midValue,maxValue);
		maxValue = sMax;
		midValue = sMid;
		minValue = sMin;
	}while(isChange && p < m);
	return 0;
}
